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In the standard perturbation theory (SPT) of self-gravitating Newtonian fluid in an expanding universe, re¬ 
currence relations for higher-order solutions are well known and play an important role both in practical ap¬ 
plications and in theoretical investigations. The recurrence relations in Lagrangian perturbation theory (LPT), 
however, have not been known for a long time. Recently, two different kinds of recurrence relations in LPT 
have been proposed in limited cases. In this paper, we generalize those methods, and most generally derive the 
recurrence relations, which are capable of including any initial condition in general models of cosmology. The 
fastest-growing modes in the general relations are identified, and simplified recurrence relations with accurate 
approximation for the time dependence are obtained. 

PACS numbers: 98.B0.-k, 98.65.-r, 04.25.Nx 


I. INTRODUCTION 

The large-scale structure (LSS) of the Universe is an important source of information about our Universe. Theoretically 
understanding physical origins and statistical properties of LSS is quite important in cosmology. The driving force for the 
evolution of LSS is gravitational instability. Initial density fluctuations are amplified by the attractive force of gravity and cause 
the origin of the present structure. The evolution of LSS is complicated for its nonlinear nature of dynamics. 

The nonlinear structure formation is complicated, and the nonlinear dynamics in general is quite hard to analytically under¬ 
stand. When the nonlinearity is considered weak, one can apply the perturbation theory. The linear theory of gravitational 
instability describes the evolution at sufficiently early stages or on sufficiently large scales of LSS. The linear theory has been 
quite successful in cosmology. Most observational signals from cosmic microwave background (CMB) radiation are understood 
by the linear theory, because the CMB radiation is emitted at a sufficiently early time, when the amplitude of density fluctuations 
are extremely small. While the amplitudes of density fluctuations in the present Universe are large on small scales, they are still 
small on large scales. 

The linear theory of density fluctuations is relatively simple to analyze, as each Fourier mode independently evolves with time. 
However, complicated couplings of mode take place when the fluctuations becomes nonlinear. Early stages of the mode coupling 
can be analytically described by higher-order perturbation theory beyond the linear approximation fit]. As recent observations 
of LSS are large enough, precise descriptions of a weakly nonlinear regime are important for cosmological analyses. In this 
respect, the higher-order perturbation theory attracts much attention these days. 

A straightforward way of describing non I inear perturbations is the standard perturbation theory (SPT), in which all the pertur¬ 
bation variables are expanded in Eulerian space [2|]. Along with recent interests in cosmological perturbation theory, theoretically 
various ways of improving the SPT have been proposed, such as the renormalized perturbation theory (RPT) 12101, effective 
field theory of large-scale structure (EFTofLSS) |5. 6|. and many others H. The SPT and its extensions fall into a category 
of Eulerian perturbation theory. An alternative to the Eulerian perturbation theory is provided by the Lagrangian perturbation 
theory (LPT), in which all the perturbation variables are expanded in Lagrangian space MU- The first-order LPT corresponds 
to the classic Zel’dovich approximation JH, which is a generalization of an exact solution in one-dimensional space jut]. Al¬ 
though the La gran gian variables are not directly observable, there exists a systematic way of predicting observable quantities 
from LPT |fl8l424il. 

As expected, there are merits and demerits in the SPT and LPT (see, e.g.. Refs. ll25Ll26tl ). and they are complementary to each 
other. While LPT cannot be extrapolated into a strongly nonlinear regime beyond the shell-crossing phase, there are several 
advantages over the SPT in a weakly nonlinear regime. In cosmological N-body simulations, it is convenient and customary to 
use the LPT to set up the initial conditions. Including the effects of redshift-space distortions into the higher-order perturbation 
theory is straightforward and natural in LPT lfl8[l27tl ' Most of the physical models of bias, such as the halo model | 28[| , peaks 
model |f29f |, excursion set peaks iU, etc. are Lagrangian bias, in which the bias relations are defined in Lagrangian space. It is 
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natural to use LPT to describe the evolution of biased objects in those models. In the formalism of integrated perturbation theory 
(iPT) lfl9U20l . l24l[3ltl . the redshift-space distortions and Lagrangian bias naturally fit into LPT. 

One of the striking advantages of SPT over LPT is that there are well-known recursive solutions of SPT for order-by-order 
expressions |f32j|. Those solutions are exact in the Einstein-de Sitter (EdS) model of cosmology, and approximate in most of 
the realistic models. The recursive solutions are practically crucial to give predictions of higher-order perturbation theory by 
numerical i nteg rations. Evaluating two- or higher-loop corrections to the power spectrum fully utilizes the recursive solutions of 

spt iUdMl- 

Until recently, the recursive solutions of LPT had been considered difficult to obtain, because the displacement fields in LPT 
are not irrotational in general, even though the Eulerian velocity fields are irrotational in the fastest-growing mode However, 
progress has made by a pioneering work 136] that describes how recursive solutions for the transverse part of displacement field 
are shown to be derived from the irrotational condition for the Eulerian velocity field. In the same reference, a hybrid procedure 
is proposed to derive the recursive solutions for the longitudinal part from the known recursive solutions of SPT, by using the 
SPT/LPT correspondence of perturbation kernels |[201 Ih1| . 

Most recently, in a somehow different context, recurrence relations of LPT are derived without resorting to SPT [37! [38!. 
Partly because the main target of the last formalism is to obtain exact solutions to the nonlinear fluid equations, the expansion 
parameter is a time parameter such as the scale factor, instead of the field values as in the usual LPT. This expansion scheme 
agrees with that of usual LPT for the fastest-growing mode in the EdS universe, but not for other modes and cosmological 
models. While the new expansion scheme of LPT based on a time parameter should be useful for solving particular problems, it 
is fair to say that recurrence relations in the usual expansion scheme of LPT are still lacking for general models of cosmology. 

The purpose of this paper is to generalize the above previous work and seek the recurrence relations within the usual framework 
of LPT. We show that they actually exist for any initial condition in general cosmology, and recursive solutions are explicitly 
derived. Besides, irrotational flows in Eulerian space do not need to be assumed to begin with, and it is explicitly shown that the 
fastest-growing modes of LPT automatically result in irrotational flows in Eulerian space. 

This paper is organized as follows. In Sec. m the framework of LPT is reviewed, and fundamental equations are derived and 
summarized in convenient forms for our purpose. In Sec. [Till the most general forms of recurrence relations of LPT are derived 
both in configuration space and in Fourier space. In Sec.UVl recurrence relations with a simple and accurate approximation of 
time dependence are presented. Conclusions are given in Sec.[V] In Appendix[A] an alternative form of recurrence relations for 
the longitudinal part is presented. In AppendixfBj recursive solutions up to seventh order in perturbations are explicitly given. 


II. THE LAGRANGIAN PERTURBATION THEORY 


For a pressureless Newtonian self-gravitating fluid, the equation of motion of the Eulerian comoving coordinates x(t) of a 
fluid element is given by Eh 


1 

x + 2 Hx = — t ), (1) 

a 1 

where a dot represents a Lagrangian derivative of time, H = a/a is the time-dependent Hubble parameter, V x = d/dx is a spatial 
derivative in Eulerian space, and <p{x, t ) is the gravitational potential. The gravitational potential is related to the density contrast 
6(x, t) = fix, t)/p - 1 through a Poisson equation, 

A x (p(x, t) = AnGpa 2 6(x, t), (2) 


where A x = V~ is the Laplacian operator. 

Independent field variables in LPT are the displacement field *F{q, t), which corresponds to the position difference of a fluid 
element between Lagrangian coordinates q and Eulerian coordinates x at a given time f, i.e., 

x(q,t) = q + •Fiqj). (3) 

As mass elements are distributed homogeneously in Lagrangian coordinates, the Eulerian density field p(x, t) is given by a 
continuity relation p(x, t)d 2 x = p<Pq, or fix, t ) = p/J(q, f), where 

J(q,t) — det(dx/dq) (4) 

is a Jacobian. Taking the divergence and rotation with respect to Eulerian coordinates, Eq. 0} is equivalent to a set of equations, 

V x -(x + 2Hx) = 4nGp^ - ij, (5) 

V.v x (x + 2Hx) = 0. 


( 6 ) 
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We define Jacobian matrix elements, 


OX' 

J IJ (q,t) = =d lJ +'F IJ (q,t), (7) 

dqj 

where the comma denotes spatial derivatives with respect to Lagrangian coordinates, i.e., Tfj = d'Fj/dqj. Eulerian spatial 
derivatives are given by 


d_ 

dxj 



( 8 ) 


in terms of Lagrangian spatial derivatives, where (7 1 represents matrix elements of the inverse Jacobian matrix. In terms of 
the Levi-Civita symbol Sijk, the determinant and inverse matrix are generally given by 


Or 1 ) 


_ l 

7 — ~^£ijk£pqrJipJjqJkri 

1 

ij ~ j k'jkpk'iqrJkqJpi - 


(9) 

( 10 ) 


The second equation is derived from the first equation and a standard formula (7 1 ), ; = 77/7, where 77 = dJ/dJij are cofactors 
of the Jacobian. The above formulas are quite useful in LPT |21 ]. 

In the following, we frequently use a differential operator 




(ii) 


where the partial derivative is taken with fixed Lagrangian coordinates, d/dt = d/dt\ q . One should notice that this is not a 
first-order differential operator. Lor any functions of time, A(t), B(t), C(t), the product rules of the above operator are given by 


T(AB) -T(A)B + AT(B) + 2AB, (12) 

T(ABC) = t(A) BC +A T(B) C + AB t(C) + 2ABC + IAbC + 2ABC, (13) 


and so forth. Substituting Eqs. 0, ® into Eqs. Q. ©, and using the above formulas, we have 

/ ~ 47rGp\ 

k'ijkk'pqrJipJjq I 7~ — “ I 7 f ;r + 87T Gp — 0, (14) 

Jij&jkpJqk?'Jqp ~ If (15) 

where a contraction identity = 5j P 5kq ~ &jq&k P is used to derive the second equation. With Eq. (0, the above set of 

equations fully describes the dynamical evolution of the displacement field V. Essentially equivalent equations are given in 
Refs. lfl3l.l2lll . although they assume the irrotational condition for the velocity field, V t x v = 0, which is not assumed here. It 
can be shown that the irrotational condition is compatible to the dynamical equation, Eq. ( IT5| >. and it is a consequence of the 
fastest-growing mode of LPT, as we will explicitly see in the following sections. 

The evolution equations above are nonlinear and it is hopeless to analytically find a general solution. When the scales of 
interest A are sufficiently larger than the typical scales of the displacement, I'/'l <sc A, it is useful to solve the nonlinear equations 
of motion by applying the perturbation theory, assuming the absolute values of displacement field are small enough. We expand 
the displacement field by a perturbation series, 


Y = Y_< ' I/(n) = ^ (1) + ^ (2) + ^ (3) + • • • , (16) 

n= 1 


where has the order of ( , / /n> )". 


III. RECURRENCE RELATIONS 

According to an identity V x (V x A) — V(V • A) - A A in the standard vector calculus, the displacement field is represented 
in a form, 


V = A~‘ [V(V • V) - V X (V X *F)] , 


( 17 ) 
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where V = d/dq is the spatial derivative in Lagrangian coordinates, and A 1 is the inverse operator of the Laplacian A = V • V. 
Specifically, the inverse Laplacian in configuration space, operating on a given function F(q), is represented by 


A ~ 1 F(q) = - 


1 

4n 



F(q') 


(18) 


From two kinds of spatial derivatives V ■ V and V x *1', one can reconstruct the displacement field by Eq. (ITTI) . Solutions of the 
Laplace’s equation, A!P = 0, should not be added to Eq. ( ITTI i. as we impose the statistical zero mean to the displacement field, 
('V) - 0. We call the first term in the right-hand side of Eq. ( 1 1 71 the longitudinal part, and the second term the transverse part. 
In Lagrangian dynamics, these two parts are generally coupled to each other. 


A. Longitudinal part 


Substituting Eq. ((TJ) into Eqs. (IT4t . we have 

(f - 4nGp) <F U = -e ijk e ipq 'F jj> (f - 2nGp) W Kq - ^s ijk e pqr ¥ Up W M [f - ^p) W Kr . 


(19) 


While the left-hand side linearly depends on the displacement field, the right-hand side consists of higher-order terms. The 
structure of this equation is the basis of the recurrence relations. Regarding the right-hand side as a source function, Eq. ( fl9] i has 
a form of inhomogeneous linear differential equation, 

(T-4nGp)g(t) = F(t). (20) 

The general solution of Eq. (l20l > is found by a standard method of ordinary differential equations as follows. The homogeneous 
equation of Eq. (l20l > is exactly the same as the linearized evolution equation of the Eulerian density contrast, 6+211 6 - 4nGpd = 0, 
and thus the two independent solutions of the homogeneous equation are the growing mode solution D+(t ) and decaying mode 
solution Z)_(0 IQ]]. For example, D+(t ) oc f 2/3 , I) (t) oc t~ l in the EdS model, and some analytic solutions are known for several 
cosmological models. The growing mode solution D + (t) is also known as the linear growth factor. The general solution of the 
inhomogeneous differential equation, Eq. (l2CTb . is given by 


g(t) — C\D+(t) + CiD~(f) + 



G(t,t')F(t’)dt', 


( 21 ) 


where C\ and C 2 are integration constants, f ln is the initial time, and 


D+(t)D-(f) - D-(t)D + (t’) 


G{t , f) = - 


( 22 ) 

D + {f)D.(f) - D.(t')D + (f) 

The last term of Eq. (|2H is a particular solution of the inhomogeneous equation. We denote this particular solution as an inverse 
of the differential operator. 


(f -4nGp) ‘ F(t) = J G(t,t')F(t')dt', 


(23) 


which is a linear operator. 

Using the notations above, a formal solution of Eq. © is given by 


V V = D + (t)A + + D (t)A -(f- 4nGpy l s ijk s ipq F up (f - 2nGp) W Kq + \e ijk s pqr F Up 'P hq if - ^-p] P k , 


(24) 


where A ± depends only on Lagrangian coordinates. The first two terms of the right-hand side give the linear solution of the 
longitudinal part. 


V • ¥ V) = D + (f)A+(q) + D-(t)A-(q). 


(25) 


The functions A ± (q) are determined by initial conditions. Substituting the perturbation series of Eq. (IT6l) into Eq. (l24l > and 
extracting a component of particular order n > 2, we have 


v • y(») = - 2 (f- 4 TrCpf' \s ijk s ipq 'vyy (t - 2ttGp) F[ m y\ 

m\+rri 2 =n 


1 

2 


2 (r - irtGpy' 


m\+m2+ni7,=n 


Sij k e pqr 'Fy' ) r: 


(m 2 ) 

iA 


r 



m(> n A 

k,r 


(26) 
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The right-hand side consists of only lower-order perturbations up to order n — 1. When the lower-order solutions ..., 
are known, the longitudinal part of 7?th order displacement field is obtained from this equation. The above Eq. (l26l > alone is 
not sufficient to determine the displacement field, because it lacks the transverse part. On the right-hand side, both longitudinal 
and transverse parts of lower-order perturbations enter, and thus both parts of displacement field are coupled to each other. We 
need a similar equation for the transverse part, which will be derived in the next subsection. There is an alternative but equivalent 
representation of Eq. (l26l >. which is described in AppendixlAl 

Because products of Levi-Civita symbol can be represented by products of the Kronecker delta, Eq. ( l26l ) can be equivalently 
expressed without Levi-Civita symbols. Specifically, the identities 



Sin S in 


dip 5 iq 6 ir 


u jp u m 
d kp 5 kq 

? jk^pqr — 

Sjp djq dj r 

dkp d k q d kr 


imply 

BijkEipqA jp B kq = TrA Tr« - Tr(AB), 

SijkSpqyAipBjqCkr = TrA Tr B TrC - TrA Tr(BC) - TrB Tr(CA) - TrC Tr (AB) + Tr (ABC) + Tr (ACB). 
Substituting these identities into Eq. (l26l . Levi-Civita symbols can be eliminated in the expression. 


(27) 


(28) 

(29) 


B. Transverse part 


Noting Eq. (HI has the form of JijXj = 0, and multiplying the inverse Jacobian matrix from the left, we have X ,■ = 0. 
Therefore Eq. (fl5l > is equivalent to 

e ijk JpjTJpk =0. (30) 

Substituting Eq. (|7} into the above equation, we have 

fVx'F = V^xf'Vf',. (31) 


An apparent solution of the homogeneous equation Tg(t) = 0 is a constant, and another solution is 

„ ,, C°° dt r°° da 

(t) “ I « 2 “ Jam c?H’ 


(32) 


which is a decaying function of time. In the EdS universe, E it) oc t l/ 3 . From the above solutions of the homogeneous equation, 
an inverse operator of T is constructed as 


J a 2 (t') [E_(t) - £_(f')] F(f)dt' = J a 2 (t') 

The formal solution of Eq. (PTT) is then given by 

V x V = B 0 + E.(t)B _ + f~ ] (VFi x TVYi), (34) 

where B^ Jq) are integration constants which depend only on Lagrangian coordinates and are divergence-free, V • By = 0. The 
first two terms of the right-hand side give the linear solution of the transverse part, 

Vxf (1) = B 0 (q) + E-(t)B-(q). (35) 

The divergence-free vectors Bo^(q) are determined by initial conditions. Substituting the perturbation series of Eq. i l 1 6| i and 
extracting a component of particular order n > 2, we have 

V x 'F {n) = Yj t~ l x . (36) 

m\+m.2=n 


I 


dt" 


a 2 {t") 


F(t')dt'. 


(33) 


The right-hand side consists of only lower-order perturbations up to order n — 1, and thus «th-order perturbations of the transverse 
part are given only by lower-order perturbations. As in the case of Eq. (l26l >. both the longitudinal and transverse parts of lower- 
order perturbations enter in the right-hand side. 

Eqs. ( l26t and d36l >. together with Eq. ( IT71 ) of each order, 

q,(n) = A -1 |^ V . .//(«)) _ V x (v X ¥' (,!) )] , 

are the closed set of general recurrence relations from which recursive solutions of LPT can be derived. 


( 37 ) 
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C. Seed values: First-order solution and initial condition 


The seed values for the recursive solutions are given by the first-order solutions, Eqs. (1251) , (l35l >. From those equations and 
Eq. ( 1371 ) with n — 1, the general solution of first order is given by 

•P W (q, t) = A- 1 {V [D + (t)A + (q) + D.(t)A-(q)] + V X [B,(q) + £_(/)*_(«)]}. (38) 

As the functions Bo-(q) are divergence free, the set of functions A ± (q), Bo (q) has 6 degrees of freedom at each Lagrangian 
point. They are completely determined by the initial condition of displacement field, f m ) and = , / /l 1 \q, f m ). 

However, physical degrees of freedom in an initial condition for a fluid element are just 4, instead of 6. They are the initial 
values of density contrast (5; n (^) = 6(q, f m ) and peculiar velocity v m {q) = v(q, fi n ) at an initial time f m . Thus the initial conditions 
for the displacement field have physically redundant degrees of freedom. 

To resolve the redundancy, we note that the initial density contrast and peculiar velocity are given by 

Sin = -V • fir, = a 1 >m, (39) 


where a m = a(t\„) is the value of the scale factor at the initial time, and we assume that the initial time is sufficiently early and 
the density contrast is well within the linear regime. It is obvious that only the functions A ± (q) and B (q) can be determined by 
d; n and f in - and the function B {] (q ) remains undetermined. Therefore, the initial conditions for the constant-transverse mode Bq 
cannot be associated with any physical quantity. This is natural, because a time-invariant rotation of displacement field is just 
a relabeling of the Lagrangian coordinates of fluid elements, without changing the density and velocity; only the time-varying 
rotation of displacement field is physically relevant to the vorticity. Thus the function B o can be considered as the “gauge” 
degrees of freedom in the initial condition, and we can freely choose this function without affecting any physical quantity 
(essentially the same argument is found in Ref. |[39] in nonexpanding background space). This function Bo has 2 degrees of 
freedom because of the divergence-free condition, and accounts for the redundancy of the initial conditions described above. The 
simplest and most natural choice is obviously Bo — 0. The arbitrariness of this gauge mode does not affect the fastest-growing 
mode of displacement field. 

The expansion of the Universe at the initial time is well described by the EdS universe, so we have I) , (f in ) cc f in 2 ^ 3 , I) (f in ) oc 
fin -1 , E-(fin) oc f in -1/3 , and Z) + (f in )/D + (f in ) = 2/3 f in , Z)_(f in )//)_(fi n ) = -l/f in , E + (t in )/E + (ti n ) = -1/3 f in . Using these relations, 
Eq. ( [39b completely determines the functions A ± {q), B (q). As a result, Eq. (l38l > can be represented by 


[3 D + (t) 

f r 2V-V in ) 

, 2 D - (t) V 


0 E-(t) V X (V X Din) 

[5 Z) + (f m ) ’ 

[ m 3 a m H m ) 

5 D4f,„) V 

I 1 ” UinHin) 



(40) 


where H m — 2/3fj n is the Hubble parameter at the initial time, and we choose Bo - 0 as described above. The Eulerian 
counterpart of Eq. (l4Qt can be found in §15 of Ref. jj]. The first term with a linear growth factor D + is the fastest-growing 
mode. In the most general case, Eq. (l40t can be used as seed values for the recurrence relations, Eqs. (f26l) . (f36l) . and 071) . 
Keeping only the fastest-growing mode is sufficiently accurate in most practical applications. Denoting 


Sl( q, t ) 


3 D+(t) L 2 V ■ v m (q) 
5 D + (f in ) inW 3 a m H m 


(41) 


and keeping only the fastest-growing mode, the Zel’dovich approximation |[T7h . = -VA -1 <5l, is recovered. Adopting the 

ZeTdovich approximation as the seed values for the recurrence relations, Eqs. (|26| >. (|36| >. 071 ) of n > 2, the resultant recursive 
solutions are also the fastest-growing mode of higher-order perturbations. 


D. Representations in Fourier space 
Applying the Fourier transform to the displacement field, 

r (Ab 

nq,t) = J —^e ikl >nk,t), 


( 42 ) 
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Eq. (l26l) . (l36l) and (l37l) are transformed to 


k ■ V (n) (k) = -i Yj f (kixk 2 )-(f-4nGp) l [v M (ki)x{t-2nGp)v (m2) (k2)\ 

m t +m 2 =n 'Jk l2 =k 

+ \ Y f [k 1 -{k 2 xk 3 )}{t-4nGp)- l {^ m 'Xk 1 ) 

mi+m2+m3=n ^123—k v 

k x P (n Xk) = i Y f (k\xk 2 )t- i [v (m 'Xk\)-t'F^ mi Xk 2 )\, 

m, -Lm ,—ii k\7 — k 


V (m2 Xk 2 ) x\t- ^j-p\ V im Xk 2 ) 


m\+m2=n 


'P (n Xk) = ^ {£ [* • P in Xk)] - k x [k x ¥ n) {k )]], 


k 2 

where we adopt notations, 


k\... n — k i + • • • + k ii , 


/ -•/ 


d i k\ 


d k n o 


SUki-n - k) ■ 


(2tt) 3 (2tt) 3 D 

The seed values for the above recurrence relations are given by the Fourier transform of Eq. G3.i.e., 


^Xk) = ¥ 


[3 D + (t) 

2 ik ■ Ojn) 

, 2 D - (t) k 

ik-V in] 

„ EJt) ik x(kx S in ) 

, 5 D + (<in)" 

[ m 3 ClinHin J 

5 D-(t m ) 

[ m a m H m j 

E— (/in) ^in^in 


(43) 

(44) 

(45) 

(46) 

(47) 


Although the time dependence is dropped from the argument of *P < - n ' > for notational simplicity, it actually does depend on the 
time variable. 

The above equations are the most general form of recurrence relations in Fourier space. When only the fastest-growing mode 
is considered, the first-order solution is given by 


V a Xk) = ^ D(t)6 0 (k ), (48) 

where D(t ) = D+(t)/D + (to ) is the linear growth factor normalized at the present time as IXto) = 1, and do(k) = 6\(k, to) is the 
linear density contrast at the present time. Considering only the fastest-growing mode, it is possible to represent the displacement 
field of each order as 


P (n Xk,t) = ± T - f L n (k u ...,k„- 

N'. J kh .„=k 


t)So(ki) ■ ■ -6o(k„), 


(49) 


where L n are time-dependent Fourier kernels. For n = 1, we have L(k\ t) = D(t)k/k 2 . 

In the following, we define and use the functions, 

S n (k u ■ ■ • ,k n ;t) = k\... n ■ L n (k\,.. .,k n ; t), (50) 

T n {ku...,k n \t) = -k 1 ... n xL n (k u ...,k n ;t). (51) 

Substituting Eq. ( l49t into Eqs. (143b —( f431 >. the recurrence relations for the fastest-growing mode of Fourier kernels are obtained 
as 

S, ,(*!,..., k, „ t) = Y " j (ki •■•mi ^ k( m i + l)..-n) • (T - 4;rGp) ' [L mi (kik mi ;t)x(T- 2 nGp) 

k^m2 i.k m i+i? • • • i kjii 0] 


...i !ra?! 

m i +mi =n 

1 2 1-mi {kmy-(m\+m2) X ^(mi+m2+l)—«)] 


2 ^ 

m\+m2+ni2=n 


4 nG 


J- d m2^km,\+\■> • ■ • 9 k mi + m2 , t) x j ^ p J L m3 (k mi + m2 +\ ,..., k n , t) 


t n (k u . 

■ •, k n ; t) — Y 


m\+m2=n 

Ln(ku • 

■ ■ > k,4 1 )= , f. 

kl-n 


(f - 4;rGp) ' |l mi (*i,..., k m] ; t ) 

„ I 

(kv ••mi X k( mi + \y.. n )T 1 1‘)}i (ki,..., k mi , t) • r T'L m2 (k mi +\, ■.., k n , t)|, 


m i !w2! 


(52) 

(53) 


( 54 ) 
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The right-hand sides of the above equations are not symmetric with respect to their wave vectors k \,..., k n , and so are the 
perturbation kernels obtained by the recurrence relations. Those wave vectors are interchangeable in Eq. < l49t . and thus only 
symmetrized kernels are physically relevant. The symmetrized kernels are obtained from unsymmetric kernels by a symmetriza- 
tion procedure 

LT (k u . ■ •, k n -1) = i Yj L ^ k P( D- ■ ■ ■ > k m\ (55) 

ptS n 

where the summation is taken for all the possible permutations S„ of the arguments. Corresponding to the relations of Eqs. G3- 
(l29l >. there are vector identities 


(k x k') ■ (L x L') 
[k ■ (k’ x k")] [L ■ {U x L")\ 


k L k L' 
k' L k' L' ’ 

k L k L' k L" 
k' L k' L' k' ■ L" 
k" L k" • L' k" ■ L" 


(56) 

(57) 


which provide an alternative expression for Eq. (l52l> . 


IV. RECURRENCE RELATIONS WITH APPROXIMATE TIME DEPENDENCE 

The recurrence relations derived above, Eqs. (j26j>, (O. ©-©, are completely general, and are applicable to any 

background cosmology. One can consider any initial conditions, and the resulting expressions contain every growing, nongrow¬ 
ing and decaying mode in general. In practical applications, however, one is interested in the fastest-growing mode. It is also 
known that the time dependence of the fastest-growing mode in higher-order perturbations is approximately given by cc D". 
This relation is exact for the EdS universe, in which Q m — 1, 0\ = 0 , D — a. There are residual time dependencies in general 
cosmology, although they are quite small for reasonable models 1IT2I lAOi-PPijl . The reason that the residual time dependencies are 
small is explained by the structure of the evolution equation in SPT 11451 46j]. A similar argument also applies in the case of LPT 
as shown below. 

Instead of the proper time t, we can use the logarithm of linear growth factor 


t = lnZ)(f), 


(58) 


as a time variable. In terms of the new variable, the operator of the type T - cmGp , which appears in the recurrence relations, 
Eqs. (l26l > and (l36l > with a = 0,4/3,2,4, is given by 


t - anGp = H 2 f 2 


rs 2 

(3 £> m \ 

±. 

3 a Q m 

8t 2 


8t 

OO | 

to 


(59) 


where / = d In Did In a is the linear growth rate and Q m = %nGpl3H 2 is the time dependent density parameter. The linear 
growth rate is approximately given by / ^ {J^ 55 for flat models nm and / ^ for Friedman models fill. If the growth rate 
is approximated by / = in Eq. ( 1591 . all the coefficients in the recurrence relations become independent of time; nth-order 
components of Eqs. ([T~9| > and (l3TT > approximately reduce to 


cP_ ]_d_ 

dr 2 2 dr 


3 

2 


V . tpM 


- 2 SijkSipq^ 

m\+m2=n 


cP_ \d_ 

dr 2 2 dr 



1 

2 


2 

m\+m2+tn-$=n 


cP_ ld_ 

8t 2 2 8t 


8 2 18 

V x V /<n) = V X 

8 2 18 

8t 2 2 8t 

/ 1 l 

m\+m, 2 =n 

8t 2 2 8t 


m(m) 
k,r ’ 


(60) 

(61) 


The linear equations with n — 1 are homogeneous, and the general solutions are given by a superposition of independent solutions 
e T = D, e“ 3r ^ 2 = I) 2 2 for the longitudinal part, and e° = 1, e T ^ 1 - I) 1/2 for the transverse part. As a consequence, decaying 
mode functions in the present approximation are replaced by £>_(/) —> e~ 3T/2 = D~ 2 ^ 2 , E (t) —> = i) l/2 . The general 

solution for the first-order displacement field is given by Eq. ( l40l ) with these replacements. In general, the differential equations 
of Eqs. < [60b and ( |6H can be recursively solved by standard methods, e.g., using the Laplace transform, including all the modes 
of time-dependence. 
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When we are interested in the fastest-growing mode, a simple logic of induction shows that the fastest-growing solutions of 
the above equations are given by , / /( " 1 oc e" T = D". This conclusion is exact in the EdS model, because Q m = f = 1. In general 
cosmology, even if approximation / - Q^ 2 is not so accurate at the present time, it is much more accurate in most of the time 
evolution. Thus the ratio , P (n ' 1 / D" in exact solutions is extremely insensitive to background cosmology, even more than what the 
approximation / - Q'J 2 would suggest Q], 

Substituting the inferred time dependence oc D" of the fastest-growing mode into Eqs. ( f60l) . ( [6Tb . the recurrence relations 
for n > 2 with approximate time dependence reduce to 


v • = -- V m (2) 

9 /_J mil 

m\+rri 2 =n 
m i <ni2 

-1 y m (3) 

6 Zj m i> 

n 

V x , P <n) = ^ 


4»q»i2 


O. m( m i)(p( m 2 ) 

t-ijkt-ipq T k q 


(2n + 3)(n - 1) 

Aim i m 2 + m 2 m-i + m 2 mi) 


^ / j —m l m 2 m 3 

m\+m2+m 3 =n 
m\<m2<m 3 


l - 


(2n + 3)(n - 1) 


s i i k s pq M m ' ) 'F (m2) 'F ( r\ 

tjK pqr , n , q h r ■> 


mi - x 


ni\ +rri2=n 
m\<m2 


( 2 ) ( 3 ) 

where M m ’ mi and are multiplicity factors defined by 


= 


11 (m 1 = m 2 ) 

12 (mi < m 2 ) 


mL 3) „ 


1 (mi = m 2 = m 3 ) 
3 (mi = m 2 < m 3 ) 
3 (mi < m 2 = m 3 ) 
6 (mi < m 2 < m 3 ) 


(62) 

(63) 


(64) 


The seed values for the above recurrence relation are given by the fastest-growing mode of Eq. (l40l) . i.e., the Zel’dovich approx¬ 
imation, 


*p (l) = -VA^l. 


(65) 


Eq. (l62l) can also be derived from Eq. (1A4I) . an alternative expression of longitudinal recurrence relations. 

The form of Eq. (f63l) is equivalent to the recurrence relations of Ref. 0 in the EdS limit which are derived from the 
irrotational condition of velocity field in Eulerian space, V t x v = 0, while our derivation does not assume the irrotational 
condition from the beginning. It is a consequence of selecting the fastest-growing mode in LPT that the Eulerian velocity field 
is irrotational. In fact, the Eulerian irrotational condition is equivalent to Vxf = V1(xV/, which is satisfied by Eq. d63l) with 
the fastest-growing mode, oc D". Therefore, unlike most treatments of LPT in the literature, there is no need for imposing 
the irrotational condition V A x v = 0 to begin with. Solutions with irrotational Eulerian flow form a subclass of general solutions 
of LPT, and the fastest-growing mode is in this subclass. Nevertheless, rotational flows are present in more general solutions 
0 - 

The form of Eq. d62l) is also equivalent to the recurrence relations of Ref. 137] in the EdS limit, which are derived from the 
Taylor series in the scale factor. Even though the last expansion scheme is not the same as the one here, they agree with each 
other in the special case of the fastest-growing mode in the EdS model. This agreement does not apply to other cosmological 
models, because the linear growth factor is not proportional to the scale factor in general. In fact, the recurrence relations in the 
ACDM model in the last expansion scheme are different from ours |[38ll . 

In Fourier space, the nth-order kernel L„ defined by Eq. ( l49l ) is proportional to D" in the present approximation. It is natural 
to separate this simple time dependence from the kernels as L n = D"L n , and the newly defined kernel L„ is extremely insensitive 
to background cosmology. Eq. (l49l) in this case is represented by 

;n» r 

= — L n (k l ,...,k n )d 0 (k l )---6 0 (k n ). 

n\ Jk h .. n =k 


(66) 
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Defining S„ = k\... n ■ L n and T„ = -k\... n x L n , the recurrence relations of Eqs. (f52l >— (f54l> reduce to 


S n {k\,...,k n ) = y ^ 


m i +m2=n 
m\<rri2 

n\ 

~6 


M m i mi 

mi \mi\ 


Amiini 


nl y 1 

"" "fi 2 j 


M, 


w l m 2 m 3 


m\+mi+m3=n 
m <ni2<tn\ 


tri[ \m 2 \ni 2 ,! 


1 - 


{k\--m\ X + ' [k>mi (ki ,..., k mi ) x L m2 (^mj+l ? • • •» ^n)j 

' (^(mi+l)--(mi+m2) ^ ^(wii+m2+l)-"«)J 


(2/2 + 3)(/i - 1) 

4(/721 ni 2 + 2222/223 + 27232221 ) 


( 2/2 + 3) (22 - 1 ) 


j^mi (^1 9 • • • 5 + l 5 • ■ ■ ? k mi + m2 ) X Lm 3 {km 1 +/M 2+19 • • • 9 *»)](, (67) 


r„(Jfci,..., *„) = (*-!)! 


m 1 +M 2 —n 
m\<m2 


2222 ~ 2721 
2721 ! 2222 * 


(^1— mi X k( mi + iy.. n ) [L mi (*i ,..., k m ^) • I. m . ( t k nl !+ ] ...., &/i)] , 


L n (k 1 ,..., k n ) — . \k\... n S n {kh, k n )) + k\... n x T n (k \,..., &»))] ■ 

k\--n 


( 68 ) 

(69) 


The seed values are given by Li(fc) = ik/k 2 , i.e., .S’ i = 1, Ti =0. Practically, the above equations should be quite useful. The 
recurrence relations above give unsymmetric kernels, and symmetric kernels 


• ■, K) sly L„(A: p(1) ,..., k m ), (70) 

peS n 

are physically relevant quantities. The vector identities of Eq. (l56l) and (f57l) provide an alternative expression for Eq. (l67l) . For 
readers’ convenience, recursive solutions up to seventh order are explicitly given in AnnendixlBl 


V. CONCLUSIONS 

In this paper, the recurrence relations in the usual framework of LPT are derived, generalizing the previously known recurrence 
relations in LPT. The newly derived relations are self-contained within the usual framework of LPT, and applicable to any 
initial condition in a general cosmological model. The most general recurrence relations are given by Eqs. (l26ll . (l36l) . (l37l) in 
configuration space, and Eqs. (f43b — d45l) in Fourier space. The resultant recursive solutions contain not only the fastest-growing 
mode but also all the other modes, with arbitrary initial conditions. In Fourier space, the perturbation kernels of the fastest- 
growing mode satisfy the recurrence relations of Eqs. (E2H54). 

The above recurrence relations are the most general. With an approximation for the time dependence, which is very accurate 
for realistic models of cosmology, we find simplified recurrence relations for the fastest-growing mode of LPT. They are given 
by Eqs. (l62l >. ( 1631 ). ( I37| i in configuration space, and Eqs. d67l i-(l69l> in Fourier space. The time dependence is explicitly solved 
in the last cases, and the corresponding recurrence relations are purely algebraic. Practically, these relations would be the most 
handy ones for future applications. Explicit recursive solutions up to seventh order are given in AppendixlBl 

Unlike most of the previous work, the irrotational condition of Eulerian velocity field, V f x v = 0, is not assumed throughout 
this work. However, the recurrence relations for the transverse part of the fastest-growing mode, Eqs. (l63l > and (l68l) . coincide 
with those derived from the irrotational condition. This means that the irrotationality of Eulerian flow is a consequence of 
selecting the fastest-growing mode in LPT, and there is no need to impose the condition from the beginning of LPT. Our results 
for the fastest-growing mode in the limit of the EdS model are fully consistent with the previously known recurrence relations. 

A straightforward application of this work is to use them for numerical evaluations of higher-order corrections to statistical 
measures of the large-scale structure, such as the power spectrum, bispectrum, trispectrum, etc. With the machinery of iPT | 20h 
and/or CLPT [2j], there is a systematic way of calculating such kinds of statistics from the higher-order LPT, including the 
effects of redshift-space distortions, nonlocal bias, and primordial non-Gaussianity Il20lI24l I3ll 148145011 . While two- or higher- 
loop nonlinear corrections have been investigated with Eulerian perturbation theory such as SPT, the applications of LPT have 
been limited to one-loop corrections (except Ref. El in which the SPT/LPT correspondence Pi is used), apparently because 
of the lack of recursive solutions in LPT. The recurrence relations derived in this paper, especially the simplest versions of 
Eqs. (l67l )—(l69l), should change the situation in this respect. 


Acknowledgments 

I thank Cornelius Rampf for his valuable comments on the manuscript. I acknowledge support from the Japanese Ministry of 
Education, Culture, Sports, Science, and Technology, Grant-in-Aid for Scientific Research (C) Program, Grants No. 24540267 
(2012) and No. 15K05074 (2015). 












11 


Appendix A: An alternative form of recurrence relations for the longitudinal part 


There is an alternative, but equivalent form of recurrence relations for the longitudinal part. Using Eq. © and the product rule 
of Eq. ( fl3l >. one sees Eq. (IT4l i is equivalent to 

— 4jlGpj (J — 1) — SijkSlmnJiljjin^kn- (Al) 

The general solution is formally given by 

J — \ — A + D+ + A O. + SijkS/mn ~ 4jlGp^ Ji[jj m Jkm (A2) 

where A ± {q) are integration constants of time. Substituting Eqs. <0, © into this equation, we have 

v T = A+D+ + A - D ~ 2 £ijk£ilmWk ' m ~ e £ijk£,mn ' F,JFj ’ m¥k ’ n + ~ 4kG ^ l £ijk£ilmWjl¥km + £ijk£,mnWilWjm¥kn \■ (A3) 

The linear solution is given by 

V- ¥ v> = A + D + +A_T»_, 

and extracting the particular order n > 2 , we have 

V-¥' (,,) =-- V sukSum^V^-l y e,*e/m„!f r( rV (m 2 V m3) 

9 / j l J K lim j,l k,m ^ / j l J K imn i,l j,m k,n 


(A4) 


m\+m2=n 


m\+m2+iri2=n 
1 


(T - 4nGp) 


2 SukSUm^^ + J] 


j . 

m\+m2=n m\+m2+mi,=n 


p .., p , m("U)|i/('"2)(p(m 3 ) 

t-ijkt-lmn Ty T- m T kn 


(A5) 


The above equation is equivalent to Eq. (l26l >. as straightforwardly confirmed by using the product rules of Eqs. (IT 2 l i and ( 1 1 3| i. 


Appendix B: Explicit solutions up to seventh order 


Explicitly writing down the perturbation kernel of each order is straightforward by applying the recurrence relations derived 
in the main text. For readers’ convenience, we manifestly show the explicit kernel functions of the fastest-growing mode with 
approximate time dependence, up to seventh order in this Appendix. The seventh-order perturbation theory is required in 
calculating, e.g., three-loop corrections to the power spectrum, etc. As in the main text, we apply notations 

S n(k 1? ■ - • s ^n) — ^1 —n ’ L n (k i, . . . , &«), (El) 

T„(k\,k „) = x L n (k u k n ), (B2) 


and results of these functions are presented below. The total kernel is given by 


Ln(ku ..., k tl ) 


k\.. 


[k\- n S n (k \,..., k n ) + k\... n x T„(k \,..., k ,,)]. 


(B3) 


from the above functions. For simplicity, the following results are unsymmetric with respect to their arguments. The sym- 
metrization procedure, Eq. (f70t should be applied to obtain symmetric kernels. 

To present the results, it is convenient to define the following functions, 


U(ki, k 2 ) 

V(kuk 2 ,k 3 ) 
W(k u k 2 ) 


l*i x k 2 \ 2 _ j k\ ■ k 2 \ 2 

k x 2 k 2 2 ~ l *1*2 J ’ 

1*1 -(k 2 xk 3 )\ 2 _ /* 1 _* 2\ 2 _ / *2 • *3 V 

*1 2 *2 2 *2 3 l *1*2 / l *2*3 / 

(* 1 X* 2 )(*i • * 2 ) 

*1 2 *2 2 


/ *3 ■ *1 \ 2 (*l -* 2 )(*2 • *3)(*3 ■ *l) 

\ *3*1 / *1 2 *2 2 *3 2 


(B4) 

(B5) 

(B6) 
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The recurrence relations enable us to deduce perturbation kernels of arbitrary order, once the first-order solution is given. As 
we are interested in the fastest-growing solution, the first-order solution is given by the Zel’dovich approximation, L\{k) = k/k 2 , 
i.e.. 


S i(k) — 1, 7"i(A:) = 0. 


(B7) 


The higher-order solutions are recursively derived from Eqs. (l67l > and (l 68 l) . The second-order solution is given by 


S 2 (k u k 2 ) = -U(k u k 2 ), T 2 (k\,k 2 ) = 0. 


(B 8 ) 


The third-order solution is given by 


S 3 (k u k 2 ,k 3 ) = ^U(k u k 23 )S 2 (k 2 ,k 3 ) - ^V(k u k 2 ,k 3 ), 
T 3 (k\,k 2 ,k 3 ) — W(k\,k 23 )S 2 (k 2 ,k 3 ). 


(B9) 

(BIO) 


Transverse parts T„ appear only from the third order (n > 3) in the fastest-growing-mode solutions. This property is a conse¬ 
quence of the fastest-growing mode with approximate time dependence, and does not apply in general solutions ll 10 . as seen 
from Eqs. (l35l > and ( [36b . The fourth-order solution is given by 


S 4 (k \,..., k 4 ) - 


28 

IT 


U(k\, k 234 )S 3 (k 2 , k 3 , £ 4 ) - W(k\, k 234 ) ■ T 3 (k 2 ,k 3 ,k 4 ) 


17 26 

+ Y^U(k n , k 34 )S 2 (k u k 2 )S 2 (k 3 , k 4 ) - —V{k\,k 2 ,k 34 )S 2 {k 3 ,k 4 ), 


T 4 (k \,..., k 4 ) = 2 


W(k u k 234 )S 3 (k 2 , k 3 , k 4 ) + —U-t-(&i x ^ 234 ) ■ T 3 (k 2 , k 3 , k 4 ) 

krk 2 34 2 


(Bll) 

(B12) 


Substituting Eqs. (IB7b — (IB 1 Ob into the above equation, explicit solutions in the EdS universe derived in Ref. (t\ I are exactly 
reproduced. 

Solutions of LPT for fifth or even higher orders are not found in the literature. The fifth-order solution from the recurrence 
relations is given by 


45 

Ss(ki,. ,.,k 5 )= — 


U(k\, k 2343 )S 4 (k 2 , ...,ks)~ W(k t , k 234 s) ■ T 4 (k 2 ,.. .,k 3 ) 


70 

+ — S 2 (k\,k 2 ) 


U(k 12 , k 343 )S 3 (k 3 , k 4 , k 3 ) - W{k{ 2 , k 343 ) ■ T 3 (k 3 , k 4 , k 3 ) 

60 j nil \ (*! x * 2 )' ^345 

- — ( V(k\,k 2 , k 345)5 3 (k 3 , k 4 , ks) h- ,, 


k\ k 2 k 343 


[(&! x k 2 ) x £ 345 ] • T 3 (k 3 , k 4 , k 3 ) 


75, 


- —V(ki,k 23 , k 43 )S 2 (k 2 , k 3 )S 2 (k 4 , k 3 ), 


T 5 (k 3 ,.. ,,k 5 ) - 3 

+ 2 S 2 (k\, k 2 ) 


W(k u k 2345 )S 4 (k 2 , ...,k 5 )+ —Ur ——^-(£1 x ^ 2345 )' T 4 (k 2 ,.... k 3 ) 


k\"k 2343 2 

W(k\ 2 , k 343 )S 3 (k 3 , k 4 , k 3 ) h— x k 343 ) • T 3 (k 3 , k 4 , k 3 ) 

k\ 2 k 34 5 


(B13) 


(B14) 
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The sixth-order solution is given by 


22 

S 6 (k i,...,* 6 ) = — 


U{k\, *23456)5 s(*2, . . ., *6) - W(*i, *23456) ' Ts(k 2 , . . . , *6) 

U(k\2, ^3456)^4(^3, • • •, *6) _ W{k\ 2 , *3456) ' T^kj, . . . , k(,) 


43 

+ —S 2(*i, *2) 


26 


5 i(k\,ki, * 3 ) [U (* 123 , *456)5 3 (^ 4 , * 5 , k(,) - 2W(*i23, * 456 ) ■ T^k^, * 5 , k^)] 

k 123 X ^456 


*123~*456~ 


• {[*123 x Ts(k\,k2, *3)] X [*456 X T^{k^ *5, *6)]} 


39 

— V(*i, *2, *3456)54(*3, . . . , *6) + 

124 
~5 
27 


(* 1 X *2) • *3456 


k\ 2 k2 2 k?,456 Z 


[(*1 X *2) X *3456] • T^ikj, k(,) 


- 2 (*2, £3) \ V{k\ , *23, *456)5 3 (*4, *5, *6) + ^ ^ ^ ^ 56 [(*1 X *23) X *456] ' r 3 (* 4 , *5, *6)] 

*1 *23“*456 


r 6 (* 1 ,...,* 6 ) = 4 

4- 55 2 (* 1 , * 2 ) 


-1 / (*12, *34, *56)52(*1, *2)52(*3, *4)52(*5, *e), 

J 

W(k u *23456)5 s(*2, • ■ ■ , *6) + ^ (*1 X *23456) ■ ^(*2, . . . , * 6 ) 


*l”*23456 2 

W(*i2, *3456)54(*3, • ■ •, *6) H-—j-—^(*12 X *3456) ' ?4(*3, . . ., *6) 

*12 *3456 


(B15) 


(B16) 
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The seventh-order solution is given by 


91 

S i(k\,. . £7) = — 


r 7 (*i 


U(k\, £234567)56(^2, . . ., £7) - W{k\ , ^234567) ' T(,(k2, . . ., £7) 

U(ki2, ^34567)^5(^3, . .. , £7) - W(k\2, £34567) ■ Ts(k-}, ..., £7) 


217 

+ -jy$2(£i,£2) 


315 

17 


U{kn2, £4567)5 3 (£t,£ 2 , £ 3 )5 4 (£ 4 , ■ • •, £7) 

- W(£ 12 3 , £4567) ■ [5 3 (£i,£ 2 , k 3 )Tj[(k4, ..., £ 7 ) - 7 , 3 (£ 1 ,£ 2 , £ 3 )54(£5,..., £7)] 

£l23 X £4567 


£l23 1 2 3 4 5 6 7 8 9 £4567" 


{[£123 x £2, £3)] x (£4567 x 74 (£ 4 , ..., £ 7 )]} 


203 . 

-yy -j V(£i, £2, £34567)54(£3, • ■ •, £7) + 


(£1 X £2) ' £ 3 4567 


£r£2 2 £34567 Z 


[(£1 x £2) x £ 34567 ] ■ Ts(k 3 , ..., £7) 


4567 


£l £23 2 £4567 2 


[(£1 x £23) x £4567] • r 4 (£ 4 , • ■ ■ , £7)] 


,£ 7 ) = 5 


8 °5 , , J w , , , „ ,, (£ix£ 23 )-£. 

yy52(£2, £3) j ^(£1, £23, £4567)54(£4, • • •, £7) + 

490 

yy y(£l, £234, £567)5 3(£2, £3, £4)53(£5, £6, £7) 

+ 25 3 (£ 2 , £ 3 , £ 4 ) (fc ' ^ * 23 ? } ' [(£1 x £ 234 ) x £ 567 ] • r 3 (£ 5 , £ 6 , £ 7 )] 

£l £234 £567 

+ *7 X£234)-£567 /(i _ j [t234 x ^ t4)] x t £ 5fi7 x ^ *,)]} 

£l-£ 23 4 £567“ ' 

yy5 2 (£l,£2)S2(£3, £4) jy(£l2, £34, £567)53(£5, £6, £7) 

+ 27 ~ 77— [(£12 x £34) x £ 567 ] ■ Ti(k 5 , £6, £7)] 

£l2"£34 £567 

W(£i, £234567)56(£2, • • ■, £7) “i-3--r(£i x £234567) • Ts(k2 ,..., £7) 


(B17) 


£1 ”£234567 2 


■ 95 2 (£l, £ 2 )[w(£t 2 , £34567)55(£3, • • •, £7) 


£17 x £ 


12 X /c 34567 
£l2 ”£'34567 2 


(£12 x £34567) • r 5 (£ 3 ,. . . , £ 7 )]. 


+ 5 


W(£i23, £4567)5 3 (£l, £2, £3)54(£4, • ■ •, £7) 
£l23 X £4567 


£l23 £4567 2 


{(£123 x £4567) • [5 3 (£i, £2, k-ffT^k^, • • ■ ,£7) - r 3 (£i,£ 2 , £3)54(£4, • • ■, £7)] 
+ [£123 x r 3 (£i,£2,£ 3 )] • [£4567 x r 4 (£ 4 , • ■ •,£7)]} ■ 


(B18) 


Continuing this kind of calculation and writing down similar expressions for n > 8 is straightforward and not difficult thanks to 
the recurrence relations. 
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